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In this paper, we introduce a dynamical Monte Carlo algorithm for spin models in which the 
number of the spins fluctuates from zero to a given number by addition and deletion of spins 
with a probabilistic rule. Such simulations are realized with a variable-system-size ensemble, a 
mixture of canonical ensembles each of which corresponds to a system with different size. The 
weight of each component of the mixture is controlled by a penalty term and systematically 
tuned in a preliminary run in a way similar to the multicanonical algorithm. In a measurement 
run, the system grows and shrinks without violating the detailed balance condition and we can 
obtain the correct canonical averages if physical quantities is measured only when its size is equal 
to the prescribed maximum size. The mixing of Markov chain is facilitated by the fast relaxation 
at small system sizes. The algorithm is implemented for the SK model of spin glass and shows 
better performance than that of a conventional heat bath algorithm. 
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In this paper, we introchice a dynamical Monte Carlo 
algorithm for spin modelsB* In the proposed algorithm, 
the number of the spins (the system size) fluctuates from 
zero to a given number by addition and deletion of spins 
with a probabilistic rule. A key feature of the algorithm 
is the use of an artificial mixture of canonical ensem- 
bles each of which corresponds to a system with differ- 
ent number of the spins. The weight of each component 
of the mixture is controlled by a penalty term for the 
system size and systematically tuned in a preliminary 
run in a. wav similar to that in multicanonical-type algo- 
rithmsJj'Erff In the measurement run, the system grows 
and shrinks without violating the detailed balance con- 
dition and the correct canonical averages are obtained 
if we measure the physical quantities only when its size 
is equal to the prescribed maximum size. Our idea has 
many common features to other algorithms based on ar- 
tificial extensions (or compositions) of the canonical en- 
sembles, say, multicanonical algorithms [J or simulated 
tempering algorithm.© The only difference is that we use 
the mixture of the canonical ensembles with different sys- 
tem sizes instead of the one with different temperatures 
in simulated tempering algorithm. 

The proposed algorithm is viewed as a variant of the 
methods based on the step-by-step construction of an 
approximate ground state of a spin system. For exam- 
ple, consider an Ising model with random couplings. If 
we add Ising spins sequentially to the system and choose 
the direction of a new spin to lower the energy in each 
step, we can construct a candidate of the ground state 
of system. At finite temperature, however, the detailed 
balance condition is hardly satisfied, if the system is as- 
sumed to only grow and never diminish. A method to 
preserve the detailed balance condition is to consider the 

* e-mail: iba@ism.ac.jp 



addition and deletion of a spin as a Metropolis-type trial 
and allowing the system freely to grow and shrink in a 
probabilistic rule. The bi-directional change of the sys- 
tem size is also useful to avoid the trapping in metastable 
states. The size of the system, however, can vary from 
zero to a desired size, only when proper values of the 
penalty for the size are given. Without a penalty for 
the system size, the system size will be stuck around an 
uncontrolled value. At this point, we introduce the idea 
of multicanonical ensemble learning proper values of 
the penalty in a preliminary run. This idea leads to the 
proposed algorithm. 

Our main aim to develop the algorithm is to accelerate 
the slow relaxation at low temperature, which is partic- 
ularly serious in Monte Carlo simulations of randomly 
frustrated system. While fast relaxation is attained in 
the high temperature components of the extended en- 
semble in simulated tempering algorithm, we expect fast 
relaxation in the "small system size" components of the 
ensemble in the present algorithm. In the latter part 
of the paper, we will show an example in which the 
proposed algorithm successfully facilitate the mixing of 
Markov chain in a spin glass problem. 

Several works with related ideas are found in the liter- 
ature. Ensembles of polymers with variable length have 
been used in the study of lattice polymers.tf Specif- 
ically, GrassbergerEP developed an algorithm in which 
the penalty to the length^af the polymer is tuned in a 
preliminary run. WildingcP explore the coexistence re- 
gion of Lennard-Jones fluid with a modification of the 
multicanonical ensemble extended in a two-dimensional 
space spanned by a pair of linear combinations of the 
number of the particles and the energy. These studies 
are regarded as natural developments of the preceding 
works on grand canonical ensembles of the correspond- 
ing systems. While the fluctuation of the system size 
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is limited in grand canonical ensembles, it can be arbi- 
trary large in the ensembles defined with tuned penalty 
terms. The situation is somewhat different in our case, 
because grandcanonical methods are, in itself,_riot popu- 
lar for spin models. On the other hand, WongEJ gives an 
interesting suggestion on variable-system-size ensembles 
in the context of statistical information processing, ft is, 
however, a part of a discussion paper and no concrete 
implementation of his method, "sequential buildup" , is 
reported.E2P 

Here we implement the algorithm for the Sherrington 
Kirkpatrick model of spin glass and show that the idea 
is worked well in this case_ .Let us consider the energy 
function of the SK model0EJ 

JV N 

Eas i }) = - 1 ^2j2 J ^ s ^ ^e{±i}, (i) 

<=1 3=1 

Zj 

The constant N is the system size and Z is the partition 
function. The coupling constants {Jij} are mutually in- 
dependent random samples from the Gaussian distribu- 
tion with a variance l/yN and a mean zero. 

In this case, it is easy to define a variable-system-size 
ensemble. Allowing the value Si — to the spin variables 
Si, the energy function eq.|l| is, in itself, viewed as the 
energy of the variable-system-size ensemble: 

N N 

E({S i }) = --Y l ^2j ij S i S j , & €{+1,-1,0}. (3) 

j=l 3=1 

The number of spins with nonzero values n = |Sj| is 
regarded as the size of a modified system. When n = N, 
i.e., all of the variables Si take the value of 1 or —1, the 
original system is recovered. 

As is mentioned earlier, we need a proper penalty func- 
tion / of the size n, if we want the system size n to fluc- 
tuate between zero and the original size N of the system. 
We set 

N N 

E{{Si}) = -^^^s-s^/w, S t G {+1,-1,0} 

i=l 3=1 

(4) 

and choose the adequate value of f(n) at each n through 
a preliminary simulation. We use a method based on a 
histogram construction, which is similar to that used in 
entropic samplingu^ and estimate the values of f(n) that 
give nearly equal frequencies of n in the interval < 
n < N . After the proper form of the penalty function / 
is determined by the iterative procedure, we perform a 
run for the measurement of the statistics. If we measure 
the statistics only when n — N , we recover the correct 
canonical averages. 

Results from an experiment with a sample (N = 
50) is shown in Fig.l. We measure the statistics 
W J J2(ij)(^i^j) 2 ' w h ere the summation runs over all 
pairs of the spins_(We measure them directly; the real 
replica techniqueEf is not used here.). The results of 
eight runs with different initial conditions and different 
random number seeds are recorded at four temperatures 



for each algorithm described below. 

In the first panel (a), the result with the proposed 
algorithm is shown. A heat bath algorithm (Gibbs sam- 
pling algorithm) is used to simulate the modified system 
defined by eq.^. At each step of the algorithm, a spin is 
randomly chosen and we substitute the state Si of the 
spin for one of the possible states 0, +1 and —1 with a 
conditional probability defined by the energy eq.^. Here, 
we define N iterations of the step as one Monte Carlo 
Step (The maximum system size N is used in the def- 
inition instead of n. Note that the value of n varies 
in the simulation.). 200000 MCS after the initial 1000 
MCS is used for tuning the penalty f(n). The values 
of f(n) are updated in every 10000 MCS in this pre- 
liminary run. Then, samples from another 200000 MCS 
are used for the calculation of the statistics. A dotted 
line in the figure indicates a corresponding result with 
exchange_Mmile XJailo (Metropolis-coupled chain) algo- 
rithmO'^llalHI'Ea'EfP Good agreement of both results 
supports the validity of the proposed algorithm. 

For the purpose of comparison, the results with con- 
ventional algorithms are shown in the second and third 
panels (b,c). The result shown in the panel (b) corre- 
spond to a conventional heat bath algorithm. The algo- 
rithm is similar to that in (a), but the original energy 
eq.|l| is used instead of eq.|] and the value of Si is not 
allowed, i.e., the system size n is fixed to N. After the 
initial 200000 MCS is discarded, samples from another 
200000 MCS are used for the calculation of statistics. 
Clearly, the variance of the results in (b) is much larger 
than that in (a). 

In the experiment in the panel (c) , we attached a sim- 
ulated annealing procedure to the simulation in (b), i.e., 
we gradually lower the temperature in a part of the ini- 
tial 200000 MCS. We start from T = 1.0 and decrease 
the temperature with the rate 5 x 10~ 6 /MCS. After we 
reach the target temperature, we keep it through the rest 
of the simulation. The addition of the initial anneal- 
ing phase will cause escape from the shallow metastable 
states around the irreverent local minima with higher 
energies. Although the variances of the results are con- 
siderably decreased in (c) compared with that in (b), 
they are still higher than that in (a). 

In the proposed algorithm, we choose a spin randomly 
when we "add" or "remove" it. Then, there is a ten- 
dency that some spins in the system are "likely to exist" 
(Si 7^ 0) and others are not. As a result, the quenched 
disorder in the system is no longer "randomly generated" 
one, when the system size n is smaller than the original 
size NeQ This does not cause any bias of the calcu- 
lated averages in our experiment because we discarded 
the samples with n ^ N from the measurement. Our ex- 
periment shows that we can get more independent sam- 
ples with the proposed algorithm than with conventional 
algorithm, even with this loss of the samples. On the 
other hand, we can design an algorithm where spins are 
added or removed with a prescribed order, e.g., Si can 
take the value only when Sj+i = 0. With this algo- 
rithm, we can use the samples with n < N as samples 
from a randomly generated system of the size n. Such 
a modification, however, will lead to longer relaxation 
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Fig. 1. (a) The proposed algorithm, (b) A standard heat bath algorithm, (c) A standard algorithm plus annealing (see text). The 
temperature T and the statistics jj(SiSj) 2 are shown in horizontal and vertical axes respectively. In each panel, the results 

of eight independent runs with different initial conditions and different random number seeds are shown (black dots). Dotted lines 
common in (a, b, c) show the values with a exchange Monte Carlo algorithm. 



time and increases a risk of the bias of the averages. 

In this paper, we proposed a Monte Carlo algorithm 
for spin models with a variable-system-size ensemble, a 
mixture of the canonical ensembles each of which cor- 
responds to a modified system with different number of 
spins. The weight of each component of the mixture is 
controlled by a penalty term for the system size and sys- 
tematically tuned in the preliminary run. The idea is 
successfully implemented for the SK model of spin glass, 
where the variable-system-size ensemble is easily realized 
with the introduction of the fictitious value of a spin 
that takes the values ±1 in the original model. We show 
that the proposed algorithm shows a better performance 
than a conventional dynamical Monte Carlo algorithm in 
our example. 

Although the success in this paper might partly rely 
on the property of the SK model, our method is fairly 
general and formally applicable to a number of models, 
say, spin models on lattices or random networks. We can 
also treat models with Potts or continuous spins. For 
the models with non-Ising spins, we should introduce a 
set of labels, which shows that the corresponding spin is 
assumed not to interact with any other spin and external 
field when the label takes a prescribed "kill" valueEiP 

At very low temperature, the proposed algorithm re- 
duced to a greedy algorithm to search ground states 



with step-by-step construction of the system and will 
be trapped in a metastable state. In that case, a two- 
dimensional extension, e.g., a simultaneous extension in 
the energy_and Ihe system size might be effective (See 
examplesBE3'E3o of successful two-dimensional exten- 
sions.). 
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